############################################################
## IRT for ordinal data: see Johnson and Albert 1999: Sec. 7.3; Sec. 5.2.3; Sec. 4.3.2 (p. 135) for core M-H/Gibbs sampler
## Ordinal data model using UN data - estimate all at once w/ vectorized data
############################################################
# source("C:\\Users\\baileyma\\Documents\\Research\\UnitedNations\\UN_AllVotesVectorized_RandomWalkPrior_Feb2013_JCR replicationMOD.R", echo=TRUE)

# load packages		## UTILITY: rm(list = ls(all = TRUE))
	#library(pscl); library(Zelig)	# pscl: inverse gamma function (rigamma) from Simon Jackman's pscl package (http://rss.acs.unt.edu/Rdoc/library/pscl/html/igamma.html)
							# zelig: oprobit function from Imai, King, Lau's Zelig package (http://cran.r-project.org/web/packages/Zelig/index.html)
## Initialize parameters
	Continuation	= 0		## 1 if use last values as starting values
	K			= 5000	## Number of MCMC iterations
	Burn			= 0000	## Number of burn in iterations - 20000 so far
	PrintSave		= 500		## Print kk and save every PrintSave iterations
	Smooth		= 0.1		## The smaller this is, the more priors matter -- 1.0 in previous run; (I've been working w/ 0.1 and 1.0)
	Thin			= 20

##
## Load and select data (Created in UN_DataProcessing_Feb2013.R); priors, storage variables
##
AllData =read.table(file = "~/Dropbox/research/RisingPowers/Voeten_replication_files/ReplicationIdealPointEstimation/AllData_Feb2013.txt")
	CSession	= AllData[,1];	Country	= AllData[,2]
	Session	= AllData[,3];	VoteID	= AllData[,5]
	yObs		= AllData[,6];	
	CName		= AllData[,7]	## DIAGNOSTIC: object.size(yObs)	# 2665552
	ObsID 	= 1:length(Country)
	ObsN		= length(ObsID)
	RClist 	= sort(unique(VoteID))
VoteStartEnd 	= read.table(file = "~/Dropbox/research/RisingPowers/Voeten_replication_files/ReplicationIdealPointEstimation/VoteStartEnd_Feb2013.txt")
IObsMat		= as.matrix(read.table(file ="~/Dropbox/research/RisingPowers/Voeten_replication_files/ReplicationIdealPointEstimation/IObsMat_Feb2013.txt"))
VoteN			= unlist(read.table(file ="~/Dropbox/research/RisingPowers/Voeten_replication_files/ReplicationIdealPointEstimation/VoteN_Feb2013.txt"))
IndN			= unlist(read.table(file ="~/Dropbox/research/RisingPowers/Voeten_replication_files/ReplicationIdealPointEstimation/IndN_Feb2013.txt"))
CountryList		= unlist(read.table(file ="~/Dropbox/research/RisingPowers/Voeten_replication_files/ReplicationIdealPointEstimation/CountryList_Feb2013.txt"))
SessionList		= unlist(read.table(file ="~/Dropbox/research/RisingPowers/Voeten_replication_files/ReplicationIdealPointEstimation/SessionList_Feb2013.txt"))
SmoothVector	= unlist(read.table(file ="~/Dropbox/research/RisingPowers/Voeten_replication_files/ReplicationIdealPointEstimation/SmoothVector_Feb2013.txt"))
GapYear		= unlist(read.table(file ="~/Dropbox/research/RisingPowers/Voeten_replication_files/ReplicationIdealPointEstimation/GapYear_Feb2013.txt"))
NN 			= length(unique(CSession))
TT 			= length(unique(VoteID))
kk			= 0				# iteration counter
Theta 		= rep(NA, NN); 		Beta 			= rep(NA, TT)
ThetaMean		= rep(0, NN);		BetaMean		= rep(NA, TT)
ThetaMeanObs 	= 1:length(ThetaMean)
VarTheta		= rep(NA, NN);		VarBeta 		= rep(NA, TT)
ThetaStore		= matrix(NA, K/Thin, NN); 	BetaStore	= matrix(NA, K/Thin, TT); 	Gamma1Store	= BetaStore; 	Gamma2Store = BetaStore
# Very big and don't need anymore: diagnostic: Accept		= matrix(0, K+Burn, TT); 	Rstore 	= matrix(0, K+Burn, TT)
## Priors ((Johnson & Albert p. 223: 	# A ~N(0, 1); epsilonIJ ~N(0, SigmaT);  PrSigmaTJ	= Inverse Gamma)
alpha 		= 0.2; lambda = 0.1	# From Johnson and Albert, p. 165 - prior parameters of inverse gamma distribution (CHECK ON THESE)
sigmaMH		= 0.1/(2 + 1)		# Johnson and Albert, p. 135 0.1/(Number of cutpoints + 1)
LagThetaPrior 	= rep(0, NN)		# For non-Random Walk change "+ LagThetaPrior[ii]/S2ThetaPrior" to  "+ ThetaPrior/S2ThetaPrior" and use: ThetaPrior = 0.0;
S2ThetaPrior 	= 1.0				# Mean and variance prior for theta
BetaPrior 		= 0.0;			S2BetaPrior 	= 1.0				# Mean and variance prior for theta
R 			= rep(NA, TT)
UniqueCountry 	= unique(Country)
MinCountryYear	= rep(NA, length(UniqueCountry))
for (cc in 1:length(UniqueCountry)) MinCountryYear[cc] = min(Session[Country==UniqueCountry[cc]])

if(Continuation == 0) { 
## STARTING VALUES: Theta, Beta, Z ("t" in Johnson and Albert, p. 166), Gamma1, Gamma2
	TotalK 				= 0
	ThetaVector				= 0*Country
	ThetaMean[CountryList==2] 	= 3.0	# US - only second session is used for priors for all these
	ThetaMean[CountryList==20] 	= 3.0	# Canada
	ThetaMean[CountryList==200] 	= 3.0	# UK
	ThetaMean[CountryList==666] 	= 3.0	# Israel
	ThetaMean[CountryList==365] 	=-2.0	# Russia 
	ThetaMean[CountryList==265] 	=-2.0	# East Germany
	## See UN_GammaStart_June2011.R for start values of gamma
		gammaMat 			= as.matrix(read.table(file ="~/Dropbox/research/RisingPowers/Voeten_replication_files/ReplicationIdealPointEstimation/gStart_Feb2013.txt"))
		Gamma1			= gammaMat[,2]
		Gamma2			= gammaMat[,3]
		CANDGamma1			= gammaMat[,2]
		CANDGamma2			= gammaMat[,3]
	USUK					= read.table(file ="~/Dropbox/research/RisingPowers/Voeten_replication_files/ReplicationIdealPointEstimation/USUK_Feb2013.txt")
	Beta 					= 0 -1* (apply(USUK, 1, mean, na.rm=T)<2) + (apply(USUK, 1, mean, na.rm=T)>2)
	Beta[is.na(Beta)==1]		= 0
	BetaVector				= rep(Beta, VoteN)
	## Start value for z (perceived trait which is t in Johnson and Albert)
		Gamma1Vector 		= rep(Gamma1, VoteN)
		Gamma2Vector 		= rep(Gamma2, VoteN)
		gLo				= -999*(yObs==1) 			+ Gamma1Vector*(yObs==2) 	+ Gamma2Vector *(yObs==3)
		gHi				= Gamma1Vector*(yObs==1) 	+ Gamma2Vector *(yObs==2) 	+ 999           *(yObs==3)
		U 				= runif(ObsN)
		Z 				= 0 + qnorm( pnorm(gLo - 0) +  U * 	 (  pnorm(gHi- 0) - pnorm( gLo- 0) )	)
		# see Gelfand,Hills, Racine-Poon and Smith (1990, 977) or Greene 3e, 179 for truncated sampling (here: theta=0 for all individuals)
		# ySim = 1: Z lower bound = -INF & upper bd = gamma1; ySim = 2: Z lower bd = gamma1 & upper bd. = gamma3; ySim = 3: Z lower bd = gamma2 & upper bd = INF
} ## END: if(Continuation == 0) { 

if(Continuation == 1) { 
	Z		= unlist(read.table(file ="~/Dropbox/research/RisingPowers/Voeten_replication_files/ReplicationIdealPointEstimation/ZSaveRW_RandWalk_Feb2013.txt"))
	Gamma1	= unlist(read.table(file ="~/Dropbox/research/RisingPowers/Voeten_replication_files/ReplicationIdealPointEstimation/Gamma1RW_RandWalk_Feb2013.txt"))
	Gamma2	= unlist(read.table(file ="~/Dropbox/research/RisingPowers/Voeten_replication_files/ReplicationIdealPointEstimation/Gamma2RW_RandWalk_Feb2013.txt"))
	gLo		= unlist(read.table(file ="~/Dropbox/research/RisingPowers/Voeten_replication_files/ReplicationIdealPointEstimation/gLoSaveRW_RandWalk_Feb2013.txt"))
	gHi		= unlist(read.table(file ="~/Dropbox/research/RisingPowers/Voeten_replication_files/ReplicationIdealPointEstimation/gHiSaveRW_RandWalk_Feb2013.txt"))
	ThetaVector	= unlist(read.table(file ="~/Dropbox/research/RisingPowers/Voeten_replication_files/ReplicationIdealPointEstimation/ThetaVectorSaveRW_RandWalk_Feb2013.txt"))
	BetaVector	= unlist(read.table(file ="~/Dropbox/research/RisingPowers/Voeten_replication_files/ReplicationIdealPointEstimation/BetaVectorSaveRW_RandWalk_Feb2013.txt"))
	BurnKPrev	= unlist(read.table(file ="~/Dropbox/research/RisingPowers/Voeten_replication_files/ReplicationIdealPointEstimation/BurnKRW_RandWalk_Feb2013.txt"))
	TotalK	= sum(BurnKPrev)
	} ## END: if(Continuation == 1) { 

##
##	Hybrid Metropolis-Hastings/Gibbs sampler
##
StartTime = date()
while(kk< K+Burn){
kk	= kk + 1;			if (kk %% 50 ==0) print(kk)

## STEP 1: SIMULATE THETAS ("Z" in Johnson and Albert, p. 166)
	## On prior, see Gelman, Carlin, Stern and Rubin p. 260, p. 254.
# RANDOM WALK PRIOR -- allow theta mean prior to equal value of theta in that country in previous period
	LagThetaPrior 					= ThetaMean[match(CountryList +(SessionList-1)/100,  CountryList +SessionList/100)] 
	LagThetaPrior[is.na(LagThetaPrior)==1 & SessionList==20]	= 	LagThetaPrior [ThetaMeanObs[is.na(LagThetaPrior)==1 & SessionList==20] -1]
	LagThetaPrior[is.na(LagThetaPrior)==1]	= 0
		# Set prior on first year to be ideal point from second year in previous MCMC iteration	
		for(cc in 1:length(UniqueCountry)) 	if(UniqueCountry[cc]!=511 & UniqueCountry[cc]!=628) LagThetaPrior[CountryList== UniqueCountry[cc] & SessionList== MinCountryYear[cc]] =  ThetaMean[CountryList==UniqueCountry[cc] & 
										  SessionList== min(SessionList[CountryList==UniqueCountry[cc] & SessionList > MinCountryYear[cc]])]
		# Zanzibar (UniqueCountry=511) and SouthSudan (UniqueCountry=628) only have 1 year of data so can't use prior from other year
		# for(cc in 1:length(UniqueCountry)) 	LagThetaPrior[CountryList==UniqueCountry[cc] & SessionList== MinCountryYear[cc]] =  ThetaMean[CountryList==UniqueCountry[cc] & SessionList== MinCountryYear[cc]+1]
		for(nn in 1:NN)				if(CountryList[nn]!=511 & UniqueCountry[cc]!=628) if(GapYear[nn]==1) LagThetaPrior[nn] = ThetaMean[CountryList==CountryList[nn] & SessionList== max(SessionList[CountryList==CountryList[nn] & SessionList< SessionList[nn]])]
			# For countries with gap in sessions, use most recent session as LagThetaPrior

for (ii in 1:NN)	{		BstarPrior 		= c(BetaVector[IObsMat[ii, 1:IndN[ii]] ], 1/(SmoothVector[ii]*Smooth))
					Bstar 		= c(BetaVector[IObsMat[ii, 1:IndN[ii]] ], 1)
					Zstar			= c(Z[IObsMat[ii, 1:IndN[ii]] ],  LagThetaPrior[CountryList==CountryList[ii] & SessionList==SessionList[ii]] )
					VarTheta[ii] 	= 1/(t(BstarPrior)%*% Bstar)
					ThetaMean[ii]	= VarTheta[ii] *(t(BstarPrior)%*%Zstar)
				#VarTheta[ii] 	= 1/(sum(BetaVector[IObsMat[ii, 1:IndN[ii]] ]^2) + 1/S2ThetaPrior)
				#VarTheta[VarTheta==Inf]=1					## For cases in which country's only vote has Beta=0 (in BetaStart)
				#ThetaMean[ii]	= VarTheta[ii] * sum(BetaVector[IObsMat[ii, 1:IndN[ii]] ] * Z[IObsMat[ii, 1:IndN[ii]] ]) + LagThetaPrior[ii]/S2ThetaPrior
			} ## END: for (ii in 1:NN)	{
Theta			= rnorm(NN, mean= ThetaMean, sd= sqrt(VarTheta))
if(sum(is.nan(Theta)) >0) stop("Theta problem")
Theta[Theta >  5] 	= 5			
Theta[Theta < -5]  	=-5
Theta 		= (Theta - mean(Theta)) / sqrt(var(Theta))			## Scale identification
for (ii in 1:NN)	ThetaVector[IObsMat[ii, 1:IndN[ii]] ] = Theta[ii]

## STEP 2: Generate candidate gammas for updating Gamma1, Gamma2 (See Johnson and Albert, p. 136)
	# see Gelfand,Hills, Racine-Poon and Smith (1990, 977) or Greene (3e, 179) for truncated sampling
	#CANDGamma1		= Gamma1 + sigmaMH * qnorm(runif(TT)*pnorm((Gamma2 - Gamma1)/sigmaMH))
	U1 			= runif(TT)    						## TT draws from uniform distribution
	CANDGamma1		= Gamma1 + sigmaMH * qnorm(U1*pnorm((Gamma2 - Gamma1)/sigmaMH)+ pnorm((-6          - Gamma1)/sigmaMH)*(1-U1))	# Truncate gamma1 at -6 on lowend
	U 			= runif(TT)    						## TT draws from uniform distribution
	#CANDGamma2 	= Gamma2 + sigmaMH * qnorm(U + pnorm((CANDGamma1 - Gamma2)/sigmaMH)*(1-U))
	CANDGamma2 		= Gamma2 + sigmaMH * qnorm(U*pnorm((6 - Gamma2)/sigmaMH) + pnorm((CANDGamma1 - Gamma2)/sigmaMH)*(1-U))	# Truncate gamma2 at 6 on high end
	CANDGamma1Vector 	= rep(CANDGamma1, VoteN)
	CANDGamma2Vector 	= rep(CANDGamma2, VoteN)
	CANDgLo		= -99              *(yObs==1)		+ CANDGamma1Vector *(yObs==2) + CANDGamma2Vector *(yObs==3)
	CANDgHi		= CANDGamma1Vector *(yObs==1) 	+ CANDGamma2Vector *(yObs==2) + 99               *(yObs==3)
	LikelihoodVector 	=  (	pnorm(CANDgHi- BetaVector* ThetaVector) - pnorm(CANDgLo- BetaVector* ThetaVector) )	/
				   (	pnorm(    gHi- BetaVector* ThetaVector) - pnorm(    gLo- BetaVector* ThetaVector) )
	LikelihoodVector[is.nan(LikelihoodVector)==1]	= 1
	for (tt in 1:TT) 		R[tt] 	=  	prod(LikelihoodVector[VoteStartEnd[tt,1]:VoteStartEnd[tt,2]] 	) * 
								(1-pnorm((CANDGamma1[tt]- Gamma2[tt])  /sigmaMH))    / (1 - pnorm((Gamma1[tt] - CANDGamma2[tt])/sigmaMH) 	)
#	Rstore[kk,] 	= R
#	Accept[kk,] 	= (runif(TT)<R) 
	AcceptTemp 		= (runif(TT)<R) 
	Gamma1 		= AcceptTemp * CANDGamma1 + (1-AcceptTemp)*Gamma1		# OLD: Accept[kk,] instead of AcceptTemp
	Gamma2 		= AcceptTemp * CANDGamma2 + (1-AcceptTemp)*Gamma2
	if(sum(is.na(Gamma1)) + sum(is.na(Gamma2))>0) stop("gamma problem")
	Gamma1Vector 	= rep(Gamma1, VoteN)
	Gamma2Vector 	= rep(Gamma2, VoteN)
	gLo		= -99          *(yObs==1) 	+ Gamma1Vector *(yObs==2) + Gamma2Vector *(yObs==3)
	gHi		= Gamma1Vector *(yObs==1) 	+ Gamma2Vector *(yObs==2) + 99           *(yObs==3)

## STEP 3: generate latent variable ("t" in Johnson and Albert, p. 166)
U 	= runif(ObsN)
Z 	= BetaVector* ThetaVector + qnorm( 	pnorm(gLo - BetaVector* ThetaVector) +  U * (pnorm(gHi - BetaVector* ThetaVector) - pnorm( gLo- BetaVector* ThetaVector) )	)
#if(max(Z)> 100) stop("Max Z >100")
Z[Z>9]	= 9; Z[Z < -9]	= -9;

## STEP 4: generate \beta_v estimates
for (tt in 1:TT)	{	VarBeta[tt] 	= 1/(	sum(ThetaVector[VoteStartEnd[tt,1]:VoteStartEnd[tt,2]]^2)	+ 1/S2BetaPrior)
				BetaMean[tt]	= VarBeta[tt] * sum(ThetaVector[VoteStartEnd[tt,1]:VoteStartEnd[tt,2]]* Z[VoteStartEnd[tt,1]:VoteStartEnd[tt,2]] )+ BetaPrior/S2BetaPrior
			} ## END: for (tt in 1:TT)	{
Beta		= rnorm(TT, mean= BetaMean, sd= sqrt(VarBeta))
if(max(Beta)> 10) stop("Max Beta >10") #FIXME: Commented to prevent error thrown 2017-10-05
BetaVector 	= rep(Beta, VoteN)

## STEP 5: Store samples
if(kk>Burn & ((kk-Burn) %% Thin == 0)) {
    ThetaStore[(kk-Burn)/Thin,  ] 	= Theta
    BetaStore[(kk-Burn)/Thin, ] 		= Beta
    Gamma1Store[(kk-Burn)/Thin, ] 	= Gamma1
    Gamma2Store[(kk-Burn)/Thin, ] 	= Gamma2
    if (kk %% PrintSave ==0){ 		write(c(Burn, kk, TotalK), 	file ="~/Dropbox/research/RisingPowers/Voeten_replication_files/5krep/BurnKRW_RandWalk_Feb2013.txt", 			ncol=1)
								write(Z, 				file ="~/Dropbox/research/RisingPowers/Voeten_replication_files/5krep/ZSaveRW_RandWalk_Feb2013.txt", 			ncol=1)
                                write(gLo, 				file ="~/Dropbox/research/RisingPowers/Voeten_replication_files/5krep/gLoSaveRW_RandWalk_Feb2013.txt", 		ncol=1)
                                write(gHi, 				file ="~/Dropbox/research/RisingPowers/Voeten_replication_files/5krep/gHiSaveRW_RandWalk_Feb2013.txt", 		ncol=1)
                                write(ThetaVector,		file ="~/Dropbox/research/RisingPowers/Voeten_replication_files/5krep/ThetaVectorSaveRW_RandWalk_Feb2013.txt", 	ncol=1)
                                write(Theta,			file ="~/Dropbox/research/RisingPowers/Voeten_replication_files/5krep/ThetaSaveRW_RandWalk_Feb2013.txt", 	ncol=1)
                                write(BetaVector, 		file ="~/Dropbox/research/RisingPowers/Voeten_replication_files/5krep/BetaVectorSaveRW_RandWalk_Feb2013.txt", 	ncol=1)
                                write(Gamma1, 			file ="~/Dropbox/research/RisingPowers/Voeten_replication_files/5krep/Gamma1RW_RandWalk_Feb2013.txt", 	ncol=1)
                                write(Gamma2, 			file ="~/Dropbox/research/RisingPowers/Voeten_replication_files/5krep/Gamma2RW_RandWalk_Feb2013.txt", 	ncol=1)
    } ## END: 		if (kk %% PrintSave ==0){
} 	## END: if(kk>Burn & ((kk-Burn) %% Thin == 0) {
} 	## END: while(kk< K){
EndTime = date();

## Save for continuation purposes
write(c(Burn, K, TotalK), 	file ="~/Dropbox/research/RisingPowers/Voeten_replication_files/5krep/BurnKRW_RandWalk_Feb2013.txt", 			ncol=1)
write(Z, 				file ="~/Dropbox/research/RisingPowers/Voeten_replication_files/5krep/ZSaveRW_RandWalk_Feb2013.txt", 			ncol=1)
write(gLo, 				file ="~/Dropbox/research/RisingPowers/Voeten_replication_files/5krep/gLoSaveRW_RandWalk_Feb2013.txt", 		ncol=1)
write(gHi, 				file ="~/Dropbox/research/RisingPowers/Voeten_replication_files/5krep/gHiSaveRW_RandWalk_Feb2013.txt", 		ncol=1)
write(ThetaVector,		file ="~/Dropbox/research/RisingPowers/Voeten_replication_files/5krep/ThetaVectorSaveRW_RandWalk_Feb2013.txt", 	ncol=1)
write(BetaVector, 		file ="~/Dropbox/research/RisingPowers/Voeten_replication_files/5krep/BetaVectorSaveRW_RandWalk_Feb2013.txt", 	ncol=1)
write(Gamma1, 			file ="~/Dropbox/research/RisingPowers/Voeten_replication_files/5krep/Gamma1RW_RandWalk_Feb2013.txt", 	ncol=1)
write(Gamma2, 			file ="~/Dropbox/research/RisingPowers/Voeten_replication_files/5krep/Gamma2RW_RandWalk_Feb2013.txt", 	ncol=1)

# Estimates
ThetaEst 		= apply(ThetaStore, 	2, 	mean)
ThetaEstDist	= t(apply(ThetaStore, 2, quantile, probs = c(0, 0.05, 0.1, 0.5, 0.9, 0.95, 1)))
BetaEst 		= apply(BetaStore, 	2, 	mean)

ThetaSummary 	= cbind(CountryList, SessionList, IndN, ThetaEst, ThetaEstDist)
write(t(ThetaSummary), 	file ="~/Dropbox/research/RisingPowers/Voeten_replication_files/5krep/ThetaSummaryRW_RandWalk_Feb2013.txt", ncol=length(ThetaSummary[1,]))
write(t(BetaEst), 	file ="~/Dropbox/research/RisingPowers/Voeten_replication_files/5krep/BetaEstRW_RandWalk_Feb2013.txt", ncol=1)
write(t(BetaStore), 	file ="~/Dropbox/research/RisingPowers/Voeten_replication_files/5krep/BetaStoreRW_RandWalk_Feb2013.txt", ncol=length(BetaStore[1,]))
write(t(ThetaStore), 	file ="~/Dropbox/research/RisingPowers/Voeten_replication_files/5krep/ThetaStoreRW_RandWalk_Feb2013.txt", ncol=length(ThetaStore[1,]))
write(t(Gamma1Store), 	file ="~/Dropbox/research/RisingPowers/Voeten_replication_files/5krep/Gamma1StoreRW_RandWalk_Feb2013.txt", ncol=length(Gamma1Store[1,]))
write(t(Gamma2Store), 	file ="~/Dropbox/research/RisingPowers/Voeten_replication_files/5krep/Gamma2StoreRW_RandWalk_Feb2013.txt", ncol=length(Gamma2Store[1,]))

print(StartTime); print(EndTime)
print(c("Number of iterations so far: ", TotalK + Burn + K))


## Create output with vote information (rcid, gamma1est, gamma2est, beta, US vote, UK etc, resolution name, votes)
BetaEst 	= read.table("~/Dropbox/research/RisingPowers/Voeten_replication_files/5krep/BetaEstRW_RandWalk_Feb2013.txt")
Gamma1Store = read.table("~/Dropbox/research/RisingPowers/Voeten_replication_files/5krep/Gamma1StoreRW_RandWalk_Feb2013.txt")
Gamma2Store = read.table("~/Dropbox/research/RisingPowers/Voeten_replication_files/5krep/Gamma2StoreRW_RandWalk_Feb2013.txt")
Gamma1Est 	= apply(Gamma1Store, 	2, 	mean); 	Gamma2Est 	= apply(Gamma2Store, 2,	mean)		## DIAGNOSTIC: cbind(Gamma2Est, rep(NA, TT), CountriesPerVote)
VoteTally	= read.table("~/Dropbox/research/RisingPowers/Voeten_replication_files/ReplicationIdealPointEstimation/VoteTally_Feb2013.txt")
VoteCode	= read.table("~/Dropbox/research/RisingPowers/Voeten_replication_files/ReplicationIdealPointEstimation/VoteCode_Feb2013.txt")
USplus	= read.table("~/Dropbox/research/RisingPowers/Voeten_replication_files/ReplicationIdealPointEstimation/USUKRussiaBrazilChinaIndiaIsrael_Feb2013.txt")

MinSession 	= rep(NA, length(unique(VoteID)))
Matches	= rep(NA, length(unique(VoteID)))
UNResName 	= rep(NA, length(unique(VoteID)))
UNResName 	= rep(NA, length(unique(VoteID)))

VoteMatrix = cbind(1:length(t(BetaEst)), VoteTally, VoteCode, BetaEst, Gamma1Est, Gamma2Est, USplus, MinSession, Matches, UNResName)

write(t(VoteMatrix ), 	file ="~/Dropbox/research/RisingPowers/Voeten_replication_files/5krep/VoteMatrix_Aug2020.txt", ncol=length(VoteMatrix[1,]))